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Abstract 

A  comparison  of  the  temperature  distributions  in  a  proton  exchange  membrane  (PEM)  fuel  cell  between  the  parallel-flow  gas  distributors  and  the 
interdigitated  gas  distributor  has  been  discussed  in  detail.  An  electrochemical-thermal  coupled  numerical  model  in  a  five-layer  membrane-electrode 
assembly  (MEA)  is  developed.  The  temperatures  for  the  reactant  fuels  as  well  as  the  carbon  fibers  in  the  porous  electrode  are  predicted  by  using 
a  CFD  technique.  The  overpotential  across  the  MEA  is  varied  to  examine  its  effect  on  the  temperature  distributions  of  the  PEM  fuel  cell.  It  is 
found  that  both  the  fuel  temperature  and  the  carbon  fiber  temperature  are  increased  with  increasing  the  total  overpotential.  In  addition,  the  fuel  and 
carbon-fiber  temperature  distributions  are  significantly  affected  by  the  flow  pattern  that  cast  on  the  gas  distributor.  Replacing  the  parallel-flow  gas 
distributor  by  the  interdigitated  gas  distributor  will  increase  the  local  maximum  temperature  inside  the  PEM  fuel  cell. 

©  2006  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

The  majority  of  heat  source  in  a  proton  exchange  membrane 
(PEM)  fuel  cell  are  the  irreversibility  of  electrochemical  reac¬ 
tions  as  well  as  the  Joule  heat.  It  rises  the  fuel  cell  temperature 
during  operation.  The  local  temperature  distribution  inside  a 
fuel  cell  has  a  strong  impact  on  the  fuel  cell  performance  since 
it  affects  the  water- vapor  distribution  by  means  of  condensa¬ 
tion.  Insufficient  cooling  may  result  in  local  hot  spots  and  thus 
dehydrate,  shrink  or  even  rupture  the  membrane.  In  addition, 
the  kinetics  of  electrochemical  reactions  directly  depends  on 
temperature.  Therefore,  from  the  viewpoint  of  promotion  in  reli¬ 
ability,  durability,  and  performance  of  PEM  fuel  cell,  an  efficient 
thermal  management  of  a  PEM  fuel  cell  becomes  crucial.  It  not 
only  removes  the  generated  heat  from  inside  the  fuel  cell  to 
the  outside  but  also  keep  the  spatial  uniformity  of  temperature 
distribution  that  avoids  local  hot  spots. 

Some  isothermal  models  have  been  developed  for  under¬ 
standing  the  adiabatic  mechanisms  inside  a  fuel  cell.  One¬ 
dimensional  models  accounting  for  diffusive  mass  transport  and 
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electrochemical  kinetics  were  developed  by  Bernardi  and  Ver- 
brugge  [1,2]  and  Springer  et  al.  [3,4].  Two-dimensional  models 
were  developed  by  Nguyen  and  White  [5],  Fuller  and  Newman 
[6],  and  Gurau  et  al.  [7].  Most  of  these  assume  some  concentra¬ 
tion  profile  of  reactant  species  along  the  channel  except  for  [7], 
which  accounts  directly  for  convective  mass  transport.  Recently, 
Shimpalee  and  Dutta  [8]  described  a  steady  state,  isothermal, 
three-dimensional,  and  single  phase  PEM  fuel  cell  model.  It 
used  a  commercial  code  to  solve  the  complete  Navier-Stokes 
equations.  Zhou  and  Liu  [9],  Um  and  Wang  [10],  and  Hwang  et 
al.  [11]  described  a  3D  model  for  PEM  fuel  cell.  Their  results 
agree  well  with  the  experimental  observations.  It  is  note  that  the 
above-modeled  results  are  based  on  the  adiabatic  conditions. 

This  paper  extends  the  author’s  previous  work  [12]  about 
the  cathodic  thermal  behaviors  to  study  the  thermal-transport 
phenomena  in  the  entire  membrane- electrode  assembly  of  a 
PEM  fuel  cell.  A  two-equation  model  is  employed  to  resolve  the 
temperature  difference  between  the  carbon  fibers  and  the  fuels 
in  the  gas  diffusion  layer.  In  addition,  conservations  of  mass, 
momentum,  species  and  charge  couple  with  the  above  thermal 
transports  with  proper  account  of  electrochemical  kinetics  and 
fluid  dynamics.  A  computational-fluid-dynamics  methodology 
is  used  to  integrate  the  above  multi-physics  transports  of  in  a 
PEM  fuel  cell.  The  possibility  of  hot  spot  inside  a  fuel  cell  is 


1204 


J.J.  Hwang,  S.J.  Liu  /  Journal  of  Power  Sources  162  (2006)  1203-1212 


Nomenclature 


Ci 


Di 

F 


j 

kc 


kf 

ks 

M 

p 

Pi,  m 


Pi,  s 

R 

S 

T 

u 

x,y 


mole  concentration  of  the  species  i  (molm  3), 


Ci  = 


{Wi/Md/Yf  0 (Oj/Mj ) 


( P/RT ) 


diffusivity  of  species  i  (m2  s-1) 

Faraday’s  constant  (96,487  Cmol-1) 
interfacial  heat  transfer  coefficient  (volumetric) 


(Wm3  K-1) 
current  density  (A  m-2) 
transfer  current  density  (Am-3) 
thermal  conductivity  of  the  solid  phase  in  the  cat¬ 
alyst  layer  (WK-1  m-1) 

thermal  conductivity  of  the  fluid  phase 
(WK-1m-1) 


thermal  conductivity  of  the  solid  phase 
(WK-1m-1) 

molecular  weight  (kg  mol-1) 
pressure  (Pa) 

possibility  of  the  electrolyte  in  the  connection  of 
the  catalyst  layer 

possibility  of  the  catalyst  in  the  connection  of  the 
catalyst  layer 

universal  gas  constant  (W  mol-1  K-1) 
source  terms  in  the  governing  equations 
temperature  (K) 
velocity  vectors  (ms-1) 
coordinate  system,  Fig.  1  (m) 


Greek  symbols 
a  symmetric  factor 

s  porosity  (gas  diffusion  layer) 

£c  porosity  of  the  catalyst  layer 

rjtot  total  overpotential  across  the  MEA  (V) 

k  permeability  (m2) 

li  viscosity  (m  s-2) 

p  density  (kg  m- 3 ) 

crm  ionic  conductivity  of  the  membrane  phase 

(£2-1  m-1) 

<rs  electric  conductivity  of  the  catalyst 

phase(£2-1  m-1) 
r  tortuosity 

vm  volume  fraction  of  the  ionic  conductor  (elec¬ 
trolyte  phase)  in  the  catalyst  layer 
vs  volume  fraction  of  the  electric  conductor  (catalyst 

phase)  in  the  catalyst  layer 

0m  potential  of  the  ionic  conductor  (electrolyte 
phase)  (V) 

0s  potential  of  the  electric  conductor  (catalyst  phase) 

(V) 

co  mass  fraction 


Subscripts 
a  anode 

c  catalyst  or  cathode 

e  energy 


eff 

effective 

f 

fluid 

HOR 

hydrogen  oxidation  reaction 

i 

species 

j 

electricity 

m 

membrane  phase,  momentum 

0 

exchange  current  density 

ORR 

oxygen  reduction  reaction 

ref 

reference 

s 

solid,  catalyst  phase 

tot 

total 

T 

transfer  current 

assessed  by  examine  the  effects  of  total  overpotential  and  gas 
distributor  geometry  on  the  temperature  distributions  inside  the 
PEM  fuel  cell.  Moreover,  the  detailed  temperature  distributions 
inside  a  PEM  fuel  cell  provide  a  comprehensive  understanding 
of  the  mechanisms  responsible  for  thermal  pathways.  It  is  help¬ 
ful  in  the  design  of  thermal  management  of  a  low-temperature 
fuel  cell. 


2.  The  model 

A  sectional  view  of  a  typical  PEM  fuel  cell  is  given  in 
Fig.  1 .  The  model  domain  is  confined  within  a  five-layer  MEA 
(membrane-electrode  assembly),  i.e.,  two  gas  diffusion  layers 
(GDLs),  two  catalyst  layers  (CLs),  and  a  proton  exchange  mem¬ 
brane  (PEM).  Each  GDL  is  in  contact  with  a  gas  distributor. 
Humidified  hydrogen  and  saturated  air  are  supplied  to  the  anodic 
gas  distributor  and  the  cathodic  gas  distributor,  respectively.  In 
the  anodic  catalyst  layer,  hydrogen  is  consumed  to  form  protons 
that  carry  the  ionic  current  to  the  cathode.  In  the  cathodic  catalyst 
layer,  the  electrochemical  reaction  not  only  consumes  the  oxy¬ 
gen  but  also  produces  the  water.  The  hydrogen  oxidation  reaction 
(HOR)  in  the  anodic  catalyst  layer  and  the  oxygen  reduction 
reaction  (ORR)  in  the  cathodic  catalyst  layer  are,  respectively, 
expressed  as  by  the  following  equations: 

1h2  — ►  H+  +  e“  (1) 

lo2  +  H+  +  e"  -*  1h20  (2) 

Both  feeds  are  regarded  as  ideal  gases  and  are  transported 
through  diffusion  and  convection.  The  electrodes  are  treated 
as  homogeneous  porous  media  with  uniform  morphological 
properties  such  as  porosity  and  permeability.  All  gases  within 
each  of  the  electrodes  exist  as  a  continuous  phase.  In  addi¬ 
tion,  water  in  the  vapor  phase  is  considered  to  simplify  the 
model. 

As  shown  in  Fig.  2,  two  kinds  of  gas  distributors  are  exam¬ 
ined,  i.e.,  the  parallel-channel  gas  distributor  (symmetric  flow) 
and  interdigitated  gas  distributor  (asymmetric  flow). 
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Fig.  1.  Sectional  view  of  the  fuel  cell  module. 


2.1.  Governing  equations 

The  ionic  and  electronic  current  balances  in  the  PEM  and 
GDLs,  based  on  the  Ohm’s  law,  are  respectively  described  as 
the  following  equations: 

V(-<TmV0m)  =  0  (3) 

V(— crsV0s)  =  0  (4) 

where  </>m  and  </>s  are  the  membrane-phase  potentials  and  the 
catalyst-phase  potentials,  respectively.  crm  and  crs  are  the  ionic 


and  electronic  conductivities  of  the  PEM  and  GDLs,  respec¬ 
tively. 

In  the  catalyst  layer  of  a  PEM  fuel  cell,  the  porous  matrix  con¬ 
tains  two  kinds  of  solid  phases,  i.e.,  ionic  conductor  (electrolyte) 
and  electronic  conductor  (catalyst).  A  potential  difference  exists 
between  the  catalyst  and  electrolyte  to  drive  the  transfer  current 
(j*t),  keeping  the  electrochemical  reaction  continuously.  The  cur¬ 
rent  passes  through  catalyst  layer  can  be  decomposed  two  parts, 
i.e.: 

i  =  is+im  (5) 


Fig.  2.  Schematic  drawings  of  the  flow  patterns  of  (a)  parallel-channel  gas  distributor  (symmetric  flow)  and  (b)  interdigitated-channel  gas  distributor  (asymmetric 
flow). 
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where  is  and  im  are  the  currents  flowing  through  the  catalyst  and 
the  electrolyte,  respectively.  Since  the  electrodes  are  electroneu¬ 
tral  everywhere,  there  is  no  charge-buildup  in  the  catalyst  layers. 
Thus,  the  charge  conservation  is: 

V  •  i  =  0  (6) 

That  is: 


be  described  as  the  following  equation: 


(14) 


V  •  is  =  -  V  •  im  (7) 

These  two  current  components  interact  through  electrochemical 
reactions.  The  electrons  are  transferred  to  the  catalyst  from  the 
electrolyte  in  the  anodic  catalyst  layer,  and  vice  versa  in  the 
cathodic  catalyst  layer.  Application  of  Ohm’s  law  to  equation 
yields  the  current  conservation: 


•  • 

Jt,HOR  —  Jo, a 


CH2 

cH2,ref 


(15) 


where  j0, c,  and  jG,a,  are  the  cathodic  and  anodic  exchange  current 
densities,  respectively. 

The  fluid  flow  in  the  porous  media  is  described  by  the  fol¬ 
lowing  equations: 


v  •  (-<xs,effV0s)  =  -  Sj  (8) 

V  •  (  crm  effV0m)  =  Sj  (9) 

where  the  sources  terms  Sj  and  —Sj  are  the  local  transfer  cur¬ 
rent  densities  corresponds  to  the  HOR  and  ORR  in  the  anode  and 
cathode,  respectively,  creating  and  consuming  protons.  They  are 
denoted  as  j’x.hor  and  j*x,oor>  respectively,  in  the  following  dis¬ 
cussion.  crSjeff  and  am?eff  are  the  effective  electronic  and  ionic 
conductivities  of  the  catalyst  and  electrolyte,  respectively.  They 
are  modeled  as: 

O^eff  —  °s(l  —  £c)  x  x  Pi, s  (10) 

^nxeff  =  crm(l  £c)  x  x  Pi, m  (11) 

where  vs  and  vm  are  the  volume  fraction  of  the  catalyst  and 
electrolyte  in  the  catalyst  layer,  respectively,  pi ,s  and  are 
the  possibilities  of  the  catalyst  and  electrolyte  in  the  connection 
of  the  catalyst  layer,  respectively  [7,11].  It  is  noted  that  only  a 
long-range  connection  of  the  same  particles  stretch  through  the 
entire  catalyst  layer  ensures  good  conductivity. 

The  present  model  takes  into  account  two  species  in  the  anode 
(H2,  H2O),  and  three  in  the  cathode  (O2,  H2O,  N2).  The  species 
transports  based  on  the  Stefan-Maxwell  multi-component  dif¬ 
fusion  are  given  by  the  following  equations: 


pu  •  Vu  =  —V p  +  V  •  (/xVu)  +  Sm  (16) 

V(pu)  =  (17) 

where  p  is  the  density,  p  the  viscosity,  u  the  velocity  vec¬ 
tor,  and  p  is  the  pressure.  The  source  term  in  the  momentum 
equations  is  based  on  the  Darcy’s  law,  representing  an  extra 
drag  force  proportional  to  fluid  viscosity  and  velocity,  and 
inversely  proportional  to  the  permeability  of  a  porous  medium, 
i.e.,  Sm  =  —(/jl/k) u,  where  k  is  the  permeability. 

As  for  the  energy  equations,  the  two-equation  model  is  used 
to  describe  the  thermal  behaviors  in  the  gas  diffusion  layer.  The 
energy  equations  for  fluid  and  solid  phases,  respectively,  are: 

(pCp)fU  •  V7f  =  V  •  (kf,effV7f)  -  5e,GDL  (18) 

0  =  V  •  (ks5effVTs)  +  5e,GDL  +  5^2  (19) 

The  source  terms  5c,gdl  and  Sq  are  the  thermal  interaction 
between  the  solid  matrices  and  the  fluids  and  the  resistive  heat¬ 
ing,  respectively.  They  are  represented  by  5e,GDL  =  —  hw  •  ( Ts  — 
7f)  and  S q  =  /<Ts>e ff,  respectively.  hv  is  the  interfacial  heat 
transfer  coefficient  (volumetric)  between  the  solid  matrices  and 
the  reactants  fluid  in  porous  medium  [14].  The  effective  thermal 
conductivities  of  both  phases  are  respectively  defined  as: 

^s,eff  —  (1  S)ks  (70) 


I  ^ 

pil  •  V&)/  —  V  •  <  p(Dj  y  ^P/,eff,  j 

{  7=1 


M 

M~j 


(12) 


The  effective  diffusivities  of  the  species  i  in  the  porous  electrode 
follows  the  Bruggemann  model,  i.e.: 

A-,eff  =  sTDi  (13) 

The  source  terms  Si  represent  the  consumption  of  the  reac¬ 
tants  during  the  electrochemical  reaction.  It  becomes  5/  = 
Jt, orr 4^02/4^  and  Si  =  }TuORMn2/2F  in  the  cathodic  and 
anodic  catalyst  layers,  respectively.  In  the  gas  diffusion  layer 
it  is  nothing.  According  to  the  Butler- Volmer  correlation  [13], 
the  relationship  among  the  local  transfer  current  density  (jx),  the 
reactant  concentrations  (c/),  and  phase  potentials  (0S,  0m)  can 


£f,eff  =  okf  (21) 

where  ks  and  kf  are  the  thermal  conductivities  of  the  solid  matrix 
and  the  reactant  fluid,  respectively. 

In  the  catalyst  layer,  the  electrochemical  reaction  occurs  at 
the  interface  of  reactant  fluid  and  catalyst.  Physically,  the  fluid 
and  solid  phases  in  the  catalyst  layer  have  the  same  temperatures, 
i.e.: 

(pCp)fU  •  V7f  =  V  •  (kc,effV7f)  +  5e,cL  +  Sn  (22) 

7f  =  rs  (23) 
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In  the  present  model,  the  radiation  heat  flux  and  the  energy 
dissipation  due  to  Joule  heating  are  neglected.  Therefore,  the 
source  term  in  the  above  equation  can  be  represented  by 
the  overpotential  heating  by  the  electrochemical  activation, 
i.e.,  Se,cL  -  jT.HOR  x  (0m  —  0S)  in  the  anodic  catalyst  layer  and 
*Se,CL  =  j’t.orr  x  (0s  —  0m)  in  the  cathodic  catalyst  layer,  respec¬ 
tively.  The  effective  thermal  conductivity  of  the  catalyst  layer  is 
determined  by  the  following  equation  [12]: 

1 

kc  eff  —  — T-  -  (24) 

(e/( 2kc  +  kf))  +  ((1  -  e)/3kc) 

where  kc  is  the  weight- averaged  conductivity  between  the  ionic 
conductor  (such  as  Nation™)  and  the  electric  conductor  (such 
as  Pt/C). 

As  for  the  PEM,  a  typical  conduction  equation  is  employed  to 
describe  the  thermal  behavior  in  the  impermeable  material,  i.e.: 

V  •  (*mvrs)  =  0  (25) 

where  km  is  the  thermal  conductivity  of  the  PEM. 

Table  1 


Porous-electrochemical  properties  of  the  present  modeled  fuel  cell 


Description 

Unit 

Value 

Exchange  current 

Anode,  jQ,a 

Am~3 

5.0  X  103 

density 

Cathode,  jQ,c 

Am  3 

1.0  X  10“3 

Symmetric  factor 

Anode,  aa 

— 

1.0 

Cathode,  ac 

— 

0.5 

Porosity 

a-GDL,  c-GDL,  s 

— 

0.5 

a-CL,  c-CL,  sc 

— 

0.5 

Permeability 

a-GDL,  c-GDL,  k 

2 

m 

1.57  X  10“12 

a-GDL,  c-GDL,  kc 

2 

m 

1.57  x  10“12 

Tortuosity 

a-GDL,  c-GDL,  r 

— 

1.5 

a-GDL,  c-GDL,  tc 

— 

1.5 

Thermal  conductivity 

GDL,  &s,eff 

Wkg-1  K-1 

1.7 

PEM,  km 

Wkg-1  K-1 

0.5 

Anode  gas,  fcf>e ff 

Wkg"1  K-1 

0.182 

Cathode  gas, 

Wkg-1  K-1 

0.051 

Electric  conductivity 

GDL  (electronic),  as 

m-1 

300 

PEM  (ionic),  crm 

14.4 

Inlet  pressure 

Anode,  pin,a 

kPa 

1.013  x  105 

Cathode,  /?in,a 

kPa 

1.013  x  105 

Mass  flow  rate 

Anode 

kgs-1 

1.19  x  10“6 

Cathode 

kgs-1 

2.98  x  10“5 

Species  mass  fraction 

Oxygen,  coq2 

— 

21.7% 

at  cathode  inlet 
(saturated  air  at 

STP) 

Water,  co h2o,c 

— 

2.1% 

Nitrogen,  <z>n9 

— 

77.2% 

Total 

— 

100% 

Species  mass  fraction 

Hydrogen,  <z>h9 

— 

76.5% 

at  anode  inlet 
(saturated  H2  at 

STP) 

Water,  &>H2o,a 

— 

23.5% 

Total 

— 

100% 

2.2.  Boundary  conditions 

The  electrochemical  and  physical  properties  used  in  the  cal¬ 
culation  are  given  in  Table  1 .  The  temperatures  for  both  feeds 
along  with  the  surfaces  of  the  gas  distributor  are  fixed  at  298  K. 
Both  outlets  of  the  module  have  an  ambient  pressure.  The  anode 
is  supplied  with  the  humidified  hydrogen  of  mass  fractions  of 
76.5%/23.5%  for  H2/H2O.  The  cathodic  side  feeds  with  the 
saturated  air  of  21.7%/2.1%/77.2%  for  O2/H2O/N2,  where  N2 
is  considered  as  an  inert  gas  and  serves  as  diluents.  Reac¬ 
tants  delivered  to  the  cathode  and  anode  are  2.98  x  10  5  and 
1.19  x  10  6 kgs  1 ,  respectively.  The  difference  of  electronic- 
conductor  potential  (0S)  between  two  contact  surfaces  between 
the  electrodes  and  gas  distributors  represents  the  total  overpoten¬ 
tial  (rjtot)  across  the  five-layer  MEA.  The  potential  at  the  contact 
surfaces  between  the  c-GDL  and  the  current  collector  is  arbi¬ 
trarily  chosen  to  be  zero,  while  the  total  overpotential  is  used 
as  boundary  condition  at  the  anodic  current  collector.  For  the 
rest  of  the  boundaries  they  have  either  insulation  or  symmetry 
conditions. 

2.3.  Numerical  schemes 

The  governing  equations  are  numerically  solved  by  the 
control-volume-based  finite  difference  method  [15].  The  dis¬ 
cretization  procedure  ensures  conservation  of  mass,  momentum, 
energy  and  concentration  over  each  control  volume.  The  upwind 
difference  scheme  is  used  to  treat  the  velocity,  while  the  central 
difference  scheme  is  employed  to  discretize  the  fluid  density, 
temperature,  and  species  concentration.  Momentum  equations 
corresponding  to  two  coordinates  are  solved,  followed  by  a 
pressure  correction  equation  that  does  the  mass  balance.  Con¬ 
centration  transport  equations  are  solved  after  the  bulk  flow 
calculation.  Velocity  control  volumes  are  staggered  with  respect 
to  the  main  control  volumes,  and  coupling  of  the  pressure  and 
velocity  fields  is  treated  via  the  SIMPLER  pressure  correction 


Fig.  3.  Polarization-curve  comparisons  between  the  numerical  results  and  pre¬ 
vious  experimental  results. 
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1 ;  in  length,  m.  (a)  r/tot  =  0.48  V  and  (b)  ^tot  =  0.66  V. 


algorithm  [15,16].  Obtaining  new  values  for  any  desired  vari¬ 
ables,  taking  into  account  the  latest  known  estimated  values  of 
the  variable  from  the  neighboring  nodes  solves  the  set  of  the  dif¬ 
ferential  equations  over  the  entire  region.  One  iteration  process 
is  complete  when,  in  line-by-line  technique,  all  lines  in  a  direc¬ 
tion  have  been  accounted  for.  Because  of  the  large  variations 
in  the  source  terms,  under-relaxation  is  necessary  for  the  depen¬ 
dent  variables  and  the  source  terms  to  achieve  convergence.  Line 
inversion  iteration  with  typical  under-relaxation  values  of  0. 1  for 
the  velocity  term  and  pressure  correction  terms  is  incorporated 
to  the  facilitated  calculation  [17].  Solutions  are  considered  to  be 
converged  at  each  test  condition  after  the  ratio  of  residual  source 
(including  mass,  momentum,  energy  and  species)  to  the  maxi¬ 
mum  flux  across  a  control  surface  becomes  below  1.03  x  10~6. 


In  the  present  work,  all  computations  are  performed  on 
503,160  ( X  by  Y)  structured,  orthogonal  meshes.  Additional 
runs  for  the  coarser  meshes,  403,120,  and  the  finer  meshes, 
603,180,  are  taken  for  a  check  of  grid  independence.  As  shown 
in  Fig.  3,  the  maximum  discrepancies  in  the  axial  velocity  and 
oxygen  concentration  profiles  between  two  grid  sizes  of  503, 160 
and  603,180  are  only  0.6%  and  0.9%,  respectively.  In  addition, 
results  indicate  a  maximum  change  of  0.6  percent  in  wall  temper¬ 
ature  distribution  between  the  solutions  of  503,160  and  603,180 
grids.  All  above  discrepancies  are  so  small  that  the  accuracy 
of  the  solutions  on  a  503,160  grids  is  deemed  satisfactory.  The 
CPU  time  ranged  from  3  to  6  h  on  a  Pentium  IV  PC  (2.8  GHz, 
2GB  RAM)  using  Windows  XP  operating  system.  Fig.  3  shows 
a  flow  chart  of  the  present  numerical  modeling. 


ritot=0.66V,  symmetric  flow 


(b)  (C)  Min:  298 


Fig.  5.  Distributions  of  solid-phase  temperatures  and  heat  fluxes  across  the  solid  matrices:  dimensions  in  temperature,  K;  in  heat  flux,  Wm  2;  in  length,  m.  (a) 
?7tot  =  0.66  V,  symmetric  flow,  (b)  r/tot  =  0.48  V,  symmetric  flow,  and  (c)  r]to t  =  0.48  V,  asymmetric  flow. 
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Fig.  6.  Distributions  of  fluid-phase  temperatures  and  the  heat  fluxes  across  the  fluids:  dimensions  in  temperature,  K;  in  heat  flux,  W  m  2;  in  length,  m.  (a)  rjtot  =  0.66  V, 
symmetric  flow,  (b)  r) tot  =  0.48  V,  symmetric  flow,  and  (c)  r) tot  =  0.48  V,  asymmetric  flow. 


3.  Results  and  discussion 


It  is  necessary  to  verify  the  present  model  prior  to  the  discus¬ 
sion  of  the  modeling  results.  To  this  aim,  the  authors  have  made 
a  comparison  of  the  I-V  curve  between  the  present  numerical 
predictions  with  the  experimental  data  [18].  As  shown  in  Fig.  3, 
the  agreement  for  the  comparison  is  good,  indicating  that  the 
present  numerical  results  are  reliable. 


3.1.  Thermal-fluid  transports 

Fig.  4(a)  and  (b)  shows  the  flow-velocity  vectors  of  rjto t  =  0-48 
and  0.66  V,  respectively,  for  the  symmetric-flow  geometries. 
Two  marks  on  each  plot  represent  the  local  maximum  and  min¬ 
imum  velocities  in  the  module.  It  is  seen  that  the  fuel  side 
(anode)  has  flow  velocities  direct  from  the  anodic  flow  chan¬ 
nel  toward  the  anodic  catalyst  layer.  The  maximum  velocities 
occur  at  the  inlet  around  the  rib  corners.  The  flow  at  about 
the  middle  of  the  rib  surfaces  (y  =  1 .0  mm)  is  nearly  stagnant. 
In  contrast,  the  velocities  on  the  oxidant  side  direct  from  the 
cathodic  catalyst  layer  toward  to  the  flow  channels.  This  is 
because  the  oxygen  reduction  reaction  (Eq.  (2))  in  the  cathodic 
catalyst  layer  not  only  consumes  the  oxygen  but  also  produces 
the  water  vapor  that  acts  as  a  mass  source  in  the  cathodic 
layer. 

3D  mappings  of  the  solid-phase  temperatures  and  the  fluid- 
phase  temperatures  are  shown  in  Figs.  5  and  6,  respectively. 
The  corresponding  projects  show  the  vectors  of  the  heat  fluxes 
across  the  solid  matrices  (Fig.  5)  and  the  fluids  (Fig.  6).  Two  point 
marks  shown  on  each  plot  indicate  the  maximum  or  minimum 
heat  fluxes  in  the  module. 

Fig.  5(a)  and  (b)  shows  the  effect  of  total  overpoten¬ 
tial  (?7tot)  on  the  solid-phase  temperature  distributions  for  the 
symmetric-flow  geometry.  It  is  seen  that  the  solid-phase  tem¬ 
perature  in  the  cathode  is  higher  than  that  in  the  anode. 
This  is  because  the  c-CL  dissipates  more  heat  by  the  elec¬ 
trochemical  reaction  than  the  a-CL.  The  maximum  tempera¬ 
ture  occurs  at  the  middle  of  the  module  cutting  across  the 
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cathodic  catalyst  layer  (y  =  1 .0  mm).  In  addition,  the  region  near 
the  rib  surface  has  a  low  solid-phase  temperature.  Near  the 
fuel  and  oxidant  entrances,  the  solid  phase  has  low  tempera¬ 
tures  due  to  the  convection  by  the  cool  inlet  fluid.  When  the 
total  overpotential  reduces  form  r/tot  =  0.66-0.48  V,  as  shown  in 
Fig.  3,  the  solid-phase  temperature  distribution  becomes  more 
even.  The  maximum  solid-phase  temperature  reduces  from  336 
to  317  K. 

Fig.  5(c)  shows  the  results  of  the  asymmetric-flow  geometry 
with  the  same  total  overpotential  as  that  in  Fig.  5(b).  Again,  the 
minimum  solid-phase  temperatures  occur  in  the  region  near  the 
rib  surfaces  (y  =  0  and  0.9  mm).  However,  the  maximum  solid- 
phase  temperature  has  been  pushed  downstream.  The  peak  of  the 
solid-phase  temperature  is  increased  from  317  to  323  K  with  the 
asymmetric-flow  geometry  instead  of  the  symmetric-flow  geom¬ 
etry.  This  is  because  the  effect  of  overpotential  heating  near  the 
cathode  outlet  becomes  more  significant  due  to  the  concentration 
polarization  by  depleting  the  oxidant. 

Attention  is  turned  to  the  heat  fluxes  across  the  solid  matri¬ 
ces.  It  is  seen  that  the  vectors  of  heat  fluxes  direct  from  the 
catalyst  layers  toward  to  the  surfaces  of  gas  distributors.  The 
gas  distributors  act  as  heat  sinks  to  absorb  the  heat  gener¬ 
ated  by  the  electrochemical  reaction  on  the  catalyst  surfaces. 
Detailed  inspection  of  these  plots  can  get  that  the  heat  fluxes 
decrease  with  decreasing  overpotential.  In  addition,  the  asym¬ 
metric  flow  (Fig.  5(c))  dissipated  more  heat  than  the  symmetric 
flow  (Fig.  5(b)).  It  is  interesting  to  note  that  the  significant  heat 
produced  by  cathodic  overpotential  not  only  is  conducted  out 
by  the  cathodic  gas  distributor  but  also  traverses  the  PEM  and 
is  removed  by  the  anodic  gas  distributors. 


Fig.  6  shows  the  fluid-phase  temperature  distributions 
in  the  GDL  and  CL.  No  data  are  displayed  in  the  PEM 
(0.4  mm <v<  0.5  mm)  because  of  its  impermeability  in  nature. 
Again,  the  fluid-phase  temperature  in  the  cathode  is  higher  than 
that  in  the  anode.  This  is  because  the  significant  heat  gener¬ 
ation  by  the  reaction  in  the  cathodic  catalyst  layer.  A  faster 
electrochemical  kinetics  of  HOR  (i.e.,  higher  exchange  current 
density,  Table  1)  in  the  a-CL  requires  a  less  overpotential  to 
drive  the  through-flow  current  in  the  fuel  cell.  Consequently,  a 
less  overpotential  heating  in  the  anodic  catalyst  layer  results  in 
a  lower  temperature.  The  fluid-phase  temperature  in  the  cathode 
increases  downstream  for  the  asymmetric-flow  geometry  due  to 
the  heat  accumulation.  The  maximum  fluid-phase  temperature 
for  the  asymmetric-flow  geometry  is  higher  than  the  symmetric- 
flow  geometry. 

The  heat  transfer  in  the  fluid  phase  is  removed  by  both  convec¬ 
tion  and  conduction.  Clearly,  the  heat  transfer  in  the  symmetric 
flows  (Fig.  6(a)  and  (b))  is  dominated  by  conduction.  That  is 
the  heat  flux  from  the  catalyst  layer  direct  toward  to  the  flow 
channels.  As  for  the  cathode  of  the  asymmetric  flow,  the  con¬ 
vection  dominates  the  heat  transfer  mechanism  due  to  the  high 
fluid  velocity. 

3.2.  Species  transports 

Fig.  7  shows  the  effect  of  total  overpotential  on  the  concen¬ 
tration  distribution  in  the  module  for  the  parallel-flow  geometry. 
The  data  shown  in  the  plots  are  normalized  by  the  corre¬ 
sponding  inlet  concentration,  i.e.,  <z>H2/^H2,in  for  the  anode  and 
<z>o2/<z>o2,in  for  the  cathode.  Comparing  these  two  figures  reveals 
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Fig.  9.  Overpotential  distributions  in  the  catalyst  layer,  rj tot  =  0.48  V,  symmetric 
flow. 
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ric  flow. 
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that  an  increase  of  total  overpotential  increases  the  fuel  and  oxi¬ 
dant  utilization  in  both  electrodes.  At  a  fixed  total  overpotential, 
the  concentration  decreases  from  the  module  entrance  to  the 
catalyst  layer.  A  local  minimum  of  the  species  concentration 
is  found  at  the  module  middle  in  the  catalyst  layer.  It  is  further 
seen  that  the  concentration  gradient  along  the  v-direction  in  the 
cathode  is  more  significant  than  that  in  the  anode. 

3.3.  Electron  transports 

Fig.  8  shows  the  potential  distributions  of  the  catalyst  phase 
0C  (Fig.  8(a))  and  membrane  phase  0m  (Fig.  8(b))  in  the  fuel 
cell.  The  anode  has  a  higher  catalyst-phase  than  the  cathode. 
It  decreases  along  the  v-direction  to  drive  the  electric  current 
overcoming  the  Ohmic  resistance.  In  the  catalyst  layer  and 
the  PEM,  the  membrane-phase  potential  has  a  relatively  uni¬ 
form  distribution  along  y-direction.  Again,  the  difference  of  the 
membrane-phase  potential  impels  the  ionic  current  through  the 
PEM. 

In  the  catalyst  layers,  the  differences  between  the  catalyst- 
phase  potential  and  the  membrane-phase  potential  represent 
the  activation  overpotentials,  i.e.,  oSL  =  (j)c  —  0m  in  the  a-CL  and 
7?a  =  0m  —  0c  in  the  c-CL.  Figs.  9  and  10  show  a  comparison 
of  the  activation  overpotential  distributions  at  several  elevations 
(y)  between  two  different  flow  patterns.  The  activation  overpo¬ 
tentials  in  both  catalyst  layers  increases  along  the  depth  of  the 
catalyst  layer.  That  is  the  largest  activation  overpotentials  occur 
at  the  both  sides  of  the  PEM.  The  activation  overpotential  in 
the  c-CL  is  significantly  higher  than  that  in  the  a-CL  for  both 
flow  patterns.  The  highest  activation  overpotential  occurs  at  the 
module  middle  (y  =  0. 1  mm)  for  the  symmetric  flow,  and  the 
downstream  station  most  for  the  asymmetric  flow. 

Fig.  1 1  shows  the  full-field  distribution  of  the  local  current 
density  in  the  module.  It  is  seen  that  the  current  directs  form 
the  surfaces  in  contact  with  the  anodic  gas  distributors  toward 
the  surfaces  in  contact  with  the  cathodic  gas  distributors.  The 
current  density  is  high  near  the  rib  surface  corner.  It  is  further 
seen  that  the  current  density  increases  with  increasing  the  total 
overpotential  inside  the  fuel  cell. 

4.  Concluding  remarks 

This  paper  has  presented  a  comprehensive  study  on  the 
thermal-transport  behaviors  in  a  PEM  fuel  cell.  A  multi-physics 


model  is  developed  accounting  for  conservations  of  the 
co-transports  of  mass/momentum/heat/species/charge.  In  the 
gas  diffusion  layer,  a  two-equation  thermal  transport  model 
is  developed  to  cope  with  the  local  thermal  non-equilibrium 
between  the  solid  matrices  and  the  fluids.  In  the  catalyst  layer,  a 
general  energy  equation  is  derived  using  the  volume-averaging 
technique,  along  with  a  local  heat  generation  resulting  from 
electrochemical  reactions.  Results  show  that  both  the  solid- 
matrix  temperature  and  the  fluid  temperature  increase  with 
increasing  the  total  overpotential  ( r/tot ).  Under  the  same  total 
overpotential,  the  maximum  temperature  is  increased  by  replac¬ 
ing  the  symmetric-flow  geometry  with  the  asymmetric-flow 
geometry. 
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